Density effects in a bulk binary Lennard-Jones system 
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Properties of local minima as a function of density are studied in a binary Lennard- 
Jones system for kinetic equipartition temperatures T = 1.0 (normal liquid), 0.5 
(supercooled liquid), and 0.4 (glass), in reduced units. The number of different 
local minima sampled, energy, pressure, normal mode angular frequencies, mean 
distance between all pairs of local minima and partial radial distribution functions 
are presented. In agreement with previous studies by Sastry [Phys. Rev. Lett. 85, 
590 (2000)] a limiting density is found at pi = 1.06 with negative pressure, below 
which the local structure of the glass and the supercooled phases are essentially the 
same, as evidenced by the partial radial distribution functions. The mean energy, 
pressure, and normal mode frequencies have values that are practically independent 
of the temperature below p\. 

PACS numbers: 61.20.Ja, 64.60.My, 61.43.Fs 

I. INTRODUCTION 

In recent years a considerable research effort has been expended to understand the com- 
plex phenomenology of supercooled liquids and glasses. Mode-coupling theory has proved 



quite successful for supercooled 
must be accounted for explicitly. 



iquids at higher temperatures, before activated dynamics 
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It was Goldstein who first related the behaviour 



of glass formers to the underlying potential energy surface (PES). \^ In this approach, the 
dynamics are separated into vibrational motion about a minimum on the PES and transi- 
tions between local minima or 'inherent structures'. [f^Q] The PES is partitioned into basins 
of attraction surrounding the local minima, where a basin of attraction is defined as a set 
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of points that lead to the same minimum along steepest-descent pathways. Recently, it has 
been possible to establish connections between the structure, dynamics and thermodynamics 
of finite systems and the PES in some detail. H 

Following Angell,^, Q, glass-forming liquids can be classified as strong or fragile. 
Fragile systems exhibit non-Arrhenius temperature dependence of transport properties such 
as the diffusion constant, usually accompanied by a significant heat capacity peak at the glass 
transition. In contrast, strong systems exhibit Arrhenius dynamics and small or negligible 
changes in the thermodynamic properties at the glass transition. 

Computer simulation based on standard molecular dynamics (MD) or Monte Carlo (MC) 
techniques, has played a si gnifi cant role in the testing and development of models for super- 



cooled liquids and glasses. |12| Interest in more direct connections to the PES has recently 
increased, with Sastry, Debenedetti and StillingerQ] characterising 'landscape-influenced' 
and 'landscape-dominated' regimes for a binary Lennard- Jones system, in agreement with 
the instantaneous normal modes picture of Donati, Sciortino and Tartaglia. LLJ Recently, a 

the dynamics. [1$] 

In the present contribution we focus on some properties of the energy landscape in a 
binary Lennard- Jones (BLJ) mixture, namely, the effect of the density on the local minima 
sampled as a function of temperature. This topic has been previously considered by Ma- 
landro et al.Q,!^ and Sastry. 21, 3] In the former studies the volume dependence of the 
local minima appearances and disappearances were examined. Similar results were found 
by HeuerQ- 

Sastry also analysed the relationship between fragility and the PES and how 
the former depends on the bulk density. One feature of the BLJ system is that the vibra- 
tional contribution to the entropy, within the harmonic approximation, decreases with the 
potential energy. Furthermore, Sastry found a limiting density, pi = 1.08 in reduced units, 
which defines a limit of stability separating spatially heterogeneous structures (below pi) 
from more homogeneous structures (above pi). p\ was also interpreted as a density limit to 
glass formation. Here, we concentrate on this limit and consider the glass and supercooled 
liquid structures in more detail as a function of density. 
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II. METHODS 



The model we used is a binary mixture of N = 256 atoms in a cubic box with periodic 
boundary conditions, containing 205 (~ 80%) A atoms and 51 (~ 20%) B atoms interacting 
according to a Lennard- Jones pair potential of the form 



V, 



a/3 



4e 



a 3 



Va/3 
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(1) 



r a p being the distance between particles a and (5. The values of the Lennard- Jones 



param- 



eters are €aa — 1-0, £bb = 0.5, Cab = 1-5, &aa = 1-0, (Tbb = 0.88 and <jab = 0.8.|24j The 
units of distance, energy, temperature, pressure, and time were taken as oaa-, e AA, ^aa/^b 
(Ub is the Boltzmann constant), caa^aa an d gaa^I '^aa) 1 ^ 2 : , with m the mass of both A and 
B atoms. The initial density was 1.2a ^ a with a fixed cutoff of 2.ha a p along with the mini- 
mum image convention. We truncated and shifted the potential with a quadratic function, 
so that the energy and its first derivative are continuous at the cutoff value. 13L l25| 
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community as it does not 
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ar dynamics time scale. 

[lol ] However, we have recently found that there are 



crystalline minima for this system, based on face-centred-cubic A atoms with trigonal or 
square prismatic B-A coordination, j^j Previous work has shown that this system exhibits 
a significant degree of non-Arrhenius behaviour at low temperatures, i.e. it is fragile in 
Angell's terminology. Q, H Q Q 

Micro canonical molecular dynamics simulations were carried out, where the classical 
equations of motion were integrated using a Verlet algorithm. We employed 10 5 equili- 
bration steps, followed by 10 6 data collection steps with a time step of 0.003 in reduced 
units. 

In order to study the local minima the instantaneous system configuration was quenched, 
every 1000th configuration, to locate a minimum. Quenching was performed using a mod- 
ified version of Nocedal's limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) 
algorithm, j^, 12, Q| Each quench was finished using a few eigenvector-following steps 4^, |3| 
to converge the root-mean-square gradient below 10~ 7 reduced units, employing full diag- 
onalisation of the analytical Hessian matrix to ensure that the stationary points have the 
correct Hessian index (the number of negative eigenvalues). In this manner we generated a 
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sample of minima at density p = 1.2. 

For each local minimum, the density of the system was changed by 0.02 in reduced units. 
To do this we changed the box length to give the density required, without rescaling the 
coordinates of the atoms in the box. Then, the potential energy was reminimized using the 
methods described above. This process was repeated in order to follow the properties of the 
minima over a wide range of density. 

We have analysed the behaviour of the local minima at densities in the range 0.6 < p < 
1.3, and for three different initial kinetic equipartition temperatures 1.0, 0.5 and 0.4. It is 
known for this system that the glass transition temperature, T g , occurs between 0.5 and 

.4.QQ,y,r 



III. RESULTS 



In this section we present the results obtained for the behaviour of the local minima at 
different temperatures and densities. At the initial density p = 1.2, the number of different 
local minima that we found from the 1000 instantaneous configurations, N min , was 1000, 
991, and 268 at T = 1.0, 0.5, and 0.4, respectively. These values reflect the increasing 
residence times in the local minima as the temperature is decreased, with a dramatic fall in 
N m in below Tg. We found that N m i n is basically independent of the density in the range of 
densities studied. 

We have also evaluated the mean potential energy for the local minima as a function 
of density for three different initial samples at T = 1.0, 0.5 and 0.4. In Fig. ^we can see 
how the lowest value of the energy occurs around p = 1.2 for each temperature. However, 
at T = 0.4 the minimum of the curve is slightly displaced to higher density, namely, to 
p = 1.22. 

A clear change of slope appears in Fig. ^ (the system is mechanically unstable) at pi = 1.06 
for each temperature, and for p < pi the mean potential energy is practically independent 
of T and increases as the density decreases. 

Another interesting property is the variation of the _ with deasity (Fig. D . In this 
case, the pressure, P, was calculated using the virial equation |46| 
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PV = (N - 1)T + 1 < ^ E r ^ ■ F ^ >> ( 2 ) 

where = N/ p is the volume of the supercell, N is the number of atoms, = r,j — and 
Fij is the force on the zth particle due to jth particle. A factor N — 1 appears because of 
momentum conservation. We used T = to calculate the pressure for the local minima. 
The lowest value of P is reached for each temperature around p\ = 1.06. This value is 



very close to Sastry's result of pi = 1.08|2l| (his supercell contained 204 A atoms and 52 B 
atoms). For p < 1, the mean pressure is practically independent of the temperature. Fig. Q 
also shows that the pressure is nearly zero at densities corresponding to the lowest values of 
the mean energy. 

The geometric mean normal mode frequency at each local minimum is 

3N-3 

v= n (^) i/(3N - 3 \ (3) 

i=l 

where the Ui are obtained by diagonalising the Hessian matrix. A mean angular frequency, 
averaged over all the local minima, can be defined as 

< lu >= 2ix < V >, (4) 

and this quantity is plotted in Fig. |21 The frequencies are important for dynamics because 
they appear in the usual transition state theory rate constant for barrier crossing. If we 
also use the information obtained in Fig. ^we can see that for p > 1.2, which corresponds 
to positive pressures (compressed regime), < uj > increases with energy, in agreement with 
Sastry. %J However, between p ~ 1.2 and pi = 1.06, which corresponds to negative pressures 
(stretched regime), we see the opposite behaviour, i.e. < uj > decreases as the density 
decreases. Then, it increases again to reach a constant value of 13 below p ~ 1.0 independent 
of T. 

We have also calculated the mean value of the distance between all pairs of local minima 
sampled: 



dij 



1 N 



\ E | r fc r i ' (5) 

where % and j indicates the ith and jth local minima. It is clear from Fig.|2]that this quantity 
increases with the initial temperature used to generate the samples and with decreasing 
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density. On the other hand, the initial MD trajectories sample the PES less extensively at 
low temperature for the same simulation time. In the regime pi = 1.06 < p < 1.3 < dij > 
does not change much, while below pi = 1.06 it increases steadily. 

The above results were checked by studying the partial radial distribution functions 
(PRDF) for all the local minima sampled, [33] 

y / N a N a \ 



N a (N a - 1) 

and 



9aa(r) = ^rr^—u ( E E S(r - M) (6) 



y IN A N B \ 

9 - 4flW = ^(?P (r - |r « l) /' (7) 

where 8{r) is the delta function. In Fig. El we show the PRDF for A particles as a function 
of the density for reoptimised local minima generated from the run at T = 1. 

The first peak is the highest in the PRDF for all densities, indicating that the AA 
correlation takes place at the first coordination shell. However, as the density is reduced 
from p = 1.3, this peak is slightly displaced to longer distance, i.e. it reflects the lower 
density, and then reaches a maximum value at density pi = 1.06. Below this density the 
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system shows the opposite trend due to the formation of voids and fractures, 
allowing the AA distances to relax back closer to their ideal value. 

Similar behaviour is exhibited in qab in Fig. El The strong attractive interaction between 
A and B particles is reflected in the first peak. However, qbb in Fig. 0] exhibits some 
new features. The weaker attractive interaction between B particles is reflected in the 
intensities of the peaks. The BB correlation takes place mainly in the first, second and 
third coordination shells, but for p = 0.6 the first coordination shell is more important. 
Furthermore, a clear splitting appears in the second peak, indicating that the local structure 
for B particles is more ordered at low density. This is because at low density, with the 
formation of voids in the structure, a dense packing appears. A similar signature is known 
. 50] 

471 ] suggested that the limiting density, in our case pi = 1.06, is a lower limit 
to glass formation. In order to elucidate this signature we compare the PRDF averaged over 
reoptimised local minima initially sampled at T = 0.5 and T = 0.4. The results in Fig. El 
reveal a clear difference between the two samples for densities above p\ = 1.06. For the lower 
temperature, corresponding to the glassy phase, the structure is more ordered. However, 



at high pressure 
Sastry et al. 
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below pi = 1.06 the PRDF are nearly the same, so it is not possible to differentiate between 
the glass and supercooled liquid structures in this way. Similar behaviour is exhibited by gBB 
in Fig. H3 However, qab, illustrated in Fig. El does not separate the two phases, indicating 
that the A-B disorder is similar in the glass and the supercooled liquid. 

The formation of a void in one particular local minimum sampled at T = 1 is shown in 
Figure [7| The void starts to form at the limiting density pi = 1.06 and then increases when 
the density is reduced further. 



IV. CONCLUSIONS 



In this paper we have considered the effect of the density on local minima of a binary 
Lennard- Jones system sampled at kinetic equipartition temperatures T = 1.0, 0.5 and 0.4. 
A limiting density is found around p\ = 1.06, below which fractures and voids begin to form. 
The latter features cause discontinuities in various properties, in agreement with Sastry's 
results. 

In the regime pi = 1.06 < p < 1.3 the main peaks in the PRDF for AA, BB, and AB 
particles are slightly displaced to longer distance relative to p — 1.3. Below pi, fractures and 
voids begin to form and the main peaks then move in the opposite direction. In the PRDF 
for BB particles at low density, a clear splitting of the second peak appears as a consequence 
of a more ordered local structure. Furthermore, below pi, the PRDF for AA, BB and AB 
particles for the samples generated at T = 0.5 and T = 0.4 are indistinguishable, so it is not 
possible to differentiate between the glass and supercooled liquid structures at low density 
in this way. 

The lowest value of the mean pressure also occurs around pi for each temperature. Below 
this density, the mean potential energy changes almost linearly with p due to a transition 
from homogeneous local minima to inhomogeneous structures with voids. 4^| The minimum 
of the mean energy potential as a function of the density corresponds to a pressure close to 
zero. 

On the other hand, at positive pressure (compressed region), we found that the mean 
normal mode frequencies increase with energy but, at negative pressure (stretched region), 
the opposite behaviour appears. At lower density, the mean frequency is practically the 
same for samples obtained at each temperature, and reaches a limiting value as the density 
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decreases. 

The initial MD trajectories sample fewer local minima at low temperature, as expected. 
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FIG. 1: Mean potential energy (top) and pressure (bottom) as a function of the density for the 
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FIG. 3: Partial radial distribution function for A particles (top) and A and B particles (bottom) 
averaged over reoptimised local minima initially sampled at T = 1. For clarity, the curves have 
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FIG. 4: Partial radial distribution function for B particles averaged over the reoptimised local 
minima initially sampled at T = 1. As the number of B particles is relatively small, qbb is 
somewhat noisy. For clarity, the curves have been displaced vertically. 
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FIG. 5: Partial radial distribution function for A particles (top) and A and B particles (bottom) 




FIG. 6: Partial radial distribution function for B particles averaged over reoptimised local minima 
from the samples generated at T = 0.5 and T = 0.4. For clarity, the curves have been displaced 
vertically. 
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FIG. 7: Perspective view of the A (gray) and B (black) atoms for one particular local minimum 
reoptimised at different densities and initially located in an MD run with T = 1 as the kinetic 
equipartition temperature. 



